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Abstract 

We study the functional characteristics of a two-gene motif consisting of a double positive feed- 
back loop and an autoregulatory negative feedback loop. The motif appears in the gene regulatory 
network controlling the functional activity of pancreatic /3-cells. The model exhibits bistability and 
hysteresis in appropriate parameter regions. The two stable steady states correspond to low (OFF 
state) and high (ON state) protein levels respectively. Using a deterministic approach, we show that 
the region of bistability increases in extent when the copy number of one of the genes is reduced from 
two to one. The negative feedback loop has the effect of reducing the size of the bistable region. Loss 
of a gene copy, brought about by mutations, hampers the normal functioning of the /3-cells giving rise 
to the genetic disorder, maturity-onset diabetes of the young (MODY). The diabetic phenotype makes 
its appearance when a sizable fraction of the /3-cells is in the OFF state. Using stochastic simulation 
techniques, we show that, on reduction of the gene copy number, there is a transition from the monos- 
table ON to the ON state in the bistable region of the parameter space. Fluctuations in the protein 
levels, arising due to the stochastic nature of gene expression, can give rise to transitions between the 
ON and OFF states. We show that as the strength of autorepression increases, the ON— >OFF state 
transitions become less probable whereas the reverse transitions are more probable. The implications 
of the results in the context of the occurrence of MODY are pointed out.. 
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1 Introduction 



Positive and negative feedback loops are frequently-occurring motifs in gene transcription regulatory 
networks and signaling pathways [1, 2]. The components of a feedback loop are genes, proteins and 
other molecules which are connected by regulatory interactions. Depending on the components and their 
interactions, feedback loops have distinct roles in diverse regulatory systems. A regulatory interaction 
is positive (negative) if an increase in the amount or activity of one component increases (decreases) 
the amount or activity of its interaction partner. A feedback loop is positive (negative) if the number 
of repressing interactions is zero or even (odd). A large number of experiments and theoretical studies 
elucidate the major functional characteristics of feedback loops with simple structure [1, 2, 3, 4, 5, 6, 
7, 8]. Positive feedback in a gene transcription regulatory network (GTRN) tends to enhance protein 
levels whereas negative feedback favours homeostasis ,i.e., maintenance of proteins at a desired level. 
The simplest feedback loop has only one component which is thus self -regulating. For such a motif in 
a GTRN, a protein promotes / represses its own production via autoactivation / autorepression of the 
expression of its gene. A positive feedback loop with two components and two regulatory interactions 
is of two types: double negative and double positive. Again, considering a GTRN, the protein products 
of the two genes in a double negative feedback loop repress each other's synthesis. The construction of 
a synthetic circuit, the genetic toggle, is based on this motif [9]. The double positive feedback loop is 
defined by two genes, the protein products of which promote each other's synthesis. There are several 
examples of two-component positive feedback loops in natural cellular networks [1, 2], a prominent 
example being the cell division cycle, the regulatory network of which contains both double positive 
and double negative feedback loops [10]. In this case, the loops control enzymatic activity. The double 
negative feedback loop, because of its more common occurrence, has been extensively studied in contrast 
to the double positive feedback loop. 

The next stage of complexity in feedback loops involves linked positive and negative feedback loops 
[2, 11, 12, 13]. The key variables in the dynamics of feedback loops are the concentrations of the com- 
ponent molecules. In the case of a GTRN, these may be the protein concentrations. In a deterministic 
description, the time evolution of the concentrations is determined by solving a set of coupled differential 
equations, the number of equations being equal to the number of variables. In reality, the biochemical 
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Figure 1 : The two-gene network model. The protein products of the GP and GS genes activate each 
other's synthesis. There is also an autorepressive loop in which the proteins of the GS gene repress their 
own synthesis. 

events associated with gene expression are probabilistic in nature and this is reflected in the presence of 
fluctuations (termed noise) around mean protein levels [14]. A stochastic description of time evolution 
is thus more appropriate. A single positive feedback loop has a tendency to amplify noise, also the time 
taken to reach the steady state protein level is longer than that in the case of an unregulated gene [1,2]. 
Interlinking of two positive feedback loops with slow and fast dynamics results in a switch with rapid ac- 
tivation and slow deactivation times and a marked resistance to noise in the upstream signaling pathway 
[11]. Addition of a single negative feedback loop leads to rapid deactivation in the absence of the signal 
which activates the switch [12]. The combination of positive and negative feedback loops may give rise 
to excitability with transient activation of protein levels. Recent experiments suggest that competence 
development in B. subtilis is achieved via excitability [15]. 

In this paper, we study the functional characteristics of a two-gene double positive feedback loop 
coupled with autorepression of the expression of one of the genes. The major motivation for studying this 
specific motif is its presence in the GTRN controlling the pancreatic /5-cell function [16]. The hormone 
insulin is a small protein that is synthesized in the /3-cells and secreted when an increase in the blood 
glucose level is sensed. Glucose metabolism releases energy needed by cells to do useful work. Insulin 
is necessary to metabolize glucose and thereby control its level in the blood. Diabetes occurs due to an 
excessive accumulation of glucose in the blood brought about by an insufficient production of or reduced 
sensitivity to insulin. The core of the /9-cell transcriptional network consists of a double positive feedback 
loop in which the transcription factors HNF — la and HNF — Aa, belonging to the nuclear hepatocyte 
family, activate each other's synthesis. There is also some evidence that HNF — Aa autorepresses its 
own synthesis [16]. Mutations in the transcription factors HNF — la and HNF — Aa give rise to a 
type of diabetes known as maturity-onset diabetes of the young (MODY) which has an early onset with 
age less than usually 25 years. There are six different forms of MODY of which MODY 1 and MODY 
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3 are caused by mutations in the genes hnf — 4a and hnf — la respectively [17]. The structure of the 
regulatory network, of which the two genes hnf — la and hnf — 4a are integral components, is not fully 
known. A partial structure of the complex network is shown in [16, 17] involving the genes hnf — la, 
hnf — 4a, shp, hnf — 1(3, hnf — 3(3, hnf — 37, hnf — 47, and pdx — 1. The genes collectively control 
the transcription of a number of important genes involved in glucose metabolism in the (3— cell. These 
include the glucose transporter 2 (Glut — 2) gene, the glucokinase gene encoding the glycolytic enzyme 
glucokinase which acts as glucose sensor and also the insulin gene. Odom et al. [18] combined chromatin 
immunoprecipitation assays with promoter microarrays to gain insight on the regulatory circuits formed 
by hnf — la and hnf — 4a. Both the proteins are found to control the activity of a large number of target 
genes in the (3— cell. This recent finding as well as earlier experiments [16] indicate that the hnf — la 
and hnf — 4a genes play a prominent role in the pancreatic (3— cell function. Mutations in the genes give 
rise to MODY resulting in the impairment of glucose-stimulated insulin secretion. Several experiments 
[16] provide clues on the possible molecular origins of MODY. The cross-regulatory interactions between 
HNF — la and HNF — 4a are switched on as pancreatic /3-cells receive the signals to differentiate. The 
double positive feedback loop has the potential for bistability, i.e., two stable steady states. The two states 
are a basal state in which the two genes have low activity and an activated state which corresponds to high 
protein levels. The states are analogous to the OFF and ON states of a switch. Normal functioning of the 
pancreatic /9-cells requires the two-gene feedback loop to be in the ON state. The circuit operation is, 
however, vulnerable to decreased gene dosage caused by mutations (in a diploid organism each gene has 
two identical copies). Genetic disorders, termed haploinsufficiency, are known to occur due to reduced 
gene dosage resulting in decreased protein levels [19, 20, 21, 22]. Gene expression noise increases the 
probability that a protein level falls below a threshold value so that the protein amount is insufficient 
for meaningful activity. The loss of vital protein functions is responsible for the occurrence of genetic 
disorders. MODY, brought about by reduced gene dosage, is thus an example of haploinsufficiency [16]. 
We construct a mathematical model to study the dynamics of the core circuit consisting of a double 
positive feedback loop coupled with autorepression of the hnf — 4a gene. We use both deterministic and 
stochastic approaches to identify the functional features of the motif and discuss their possible relevance 
in the occurrence of MODY 
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Figure 2: The reaction kinetic scheme of the two-gene model. The meanings of the symbols are explained 
in the text. 

2 Deterministic Approach 

The circuit diagram of the motif to be studied is shown in the figure 1 . GP and GS represent the genes 
hnf — la and hnf — Act respectively. The arrow sign denotes activation by the appropriate protein 
product and the hammerhead sign denotes repression. The chemical kinetic schemes corresponding to 
the expression of genes GP and GS are shown in figures 2(a) and 2(b). The protein products of GP 
and GS are denoted by P and S . We assume that the regulation of gene expression is mediated by the 
protein dimers S2 and P2, K P and K s being the binding constants of dimerization. For each gene, there 
are two rates of protein synthesis: a basal rate ( rate constants J P0 and J S o) and an activated rate ( rate 
constants J P and J5). In the second case, protein synthesis occurs in the activated state of the gene (GP* 
and GS*) attained via the binding of protein dimers S2 and P2 to the genes GP and GS respectively. 
The associated binding constants are K PP and K S s ■ The rate constants for protein degradation are 7p 
and 7 S with <fi denoting the degradation product. Dimer degradation is not taken into account as its rate 
is few-fold lower than the degradation rate of protein monomers. For the gene GS, there is an extra 
biochemical event representing autorepression. The dimers S2 and P2 bind the promoter region of the 
gene GS competitively, i.e., the binding of one type of dimer excludes the binding of the other type. 
When the dimer S2 binds at GS, there is complete repression. The binding constant is denoted by K R . 

The protein concentrations S and P are the dynamical variables in the system. The time scale of 
binding events, in general, is much faster than that of protein synthesis and degradation. The bound 
complexes thus reach the steady state at an earlier time point. Taking this into account, the differential 
rate equations describing the time evolution of the protein concentrations S and P are: 

7 Q 

— = JsGS* + J S oGS- ls S (1) 
at 

dP 

— = J P GP* + J P0 GP- lP P (2) 
at 
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with 



as- = ^» = r = TT1 ^ ,3, 



, 2 



GP* = VJV/ 2 , AT 2 = : (4) 



There are two conservation equations for the total concentrations n s and n P of the genes G5 1 and 
GP. 

ns = GS + GS* + (5) 

np = GP + GP* (6) 
After an appropriate change in variables 

S P 

u= , v = ——. — , r = 7 5 t (7) 

J so lis Jpo/lP 

the differential rate equations (1) and (2) are transformed into 



du 1+ripv 2 

T = ns TT< 2^ i a 2 ~ u ^ 

dr (1 + jj,u z ) + p v l 



dv_^l+£au 2 ^ ^ 
dr 1 + a u 2 

The different parameters are given by 

V = ^, i = ^ »=(^) 2 K S K R , a=(^) 2 K PP K s , (3 = f^Y K ss K P (10) 
Jso Jpo \ Is J V 75 / V 7p / 

The variable r is dimensionless whereas the variables u and v have the dimensions of concentration 
expressed in units of [nm]. The parameters rj and £ are dimensionless while the parameters /i, a and (3 
are expressed in units of The dimensions of n s and n p are in units of [nm] with one gene copy 
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Figure 3: u versus 77 curves showing bistability and hysteresis. The solid (dashed) lines represent stable 
(unstable) steady states for gene copy numbers (a) rip = 2, n s = 2 and (b) np — 2, ns — 1 . The 
parameter //, a measure of the autorepression strength, is zero. 

corresponding to approximately 1 [nm] . From now on, the units will not be explicitly mentioned. Table 1 
displays all the parameters and rescaled parameters of the two-gene model as well as their meanings and 
defining formulae. 

We use the software package XPPAUT [23] to probe the dynamics of the double positive feedback 
loop and the effect of autorepression of the S proteins on the dynamics . We focus on how the steady 
state value of u (rescaled concentration of S proteins) changes as a function of the different parameters 
in equations (8) and (9). In the steady state, the rates of change ^ and ^ are zero. Figure 3(a) shows a 
plot of u versus 77 when the autorepression strength given by ix is zero. The other parameters have values 
n s = n P = 2, f = 30.0, J so = J P0 = 2.0, 75 = 7 P = 1.0 and a = f3 = 0.002857. The plot shows 
that a region of bistability separates two region of mono stability. The two stable states in the bistable 
region correspond to low and and high values of u . In this region and at a specific value of 77, the choice 
between the stable steady states is history-dependent, i.e., depends on initial conditions [24]. If the value 
of 77 is initially low, the system ends up in the low u state. As 77 increases, the system enters the region 
of bistability but continues to be in the low expression state till a bifurcation point is reached. At this 
point, a discontinuous jump to the high u state occurs and the system becomes monostable. Bistability is 
accompanied by hysteresis , i.e., the value of rj at which the switch from the low to the high expression 
state occurs is greater than the value of rj ( the lower bifurcation point ) at which the reverse transition 
takes place. The two stable branches are separated by a branch of unstable steady states (dash-dotted 
line) which are not experimentally accessible. There are now several known systems in which bistability 
and hysteresis have been observed experimentally [9, 13, 25, 26, 27, 28, 29]. Figure 3(b) shows the plot 
of u versus rj for the same parameter values as in figure 3(a) except that the copy number of the GS gene 
is reduced from two to one, i.e., n s has the value I. A comparison of figures 3(a) and (b) shows that 
with reduced copy number the extent of the region of bistability in considerably increased. The same 
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Parameter/Rescaled Parameter 

Jpo, Jso 
Jp, Js 
Kp,K s 

Kpp, K S s 

1p, Is 

K R 
V 

a 

P 

T, M, N 

Table 1: Parameters, rescaled parameters, their meanings and defining formulae 

conclusion is reached when the steady state values of u are plotted versus the parameter (3. The region of 
bistability is lower in extent when the parameter a measure of the autorepression strength, is increased 
from zero. The value of \x is changed by modifying the value of K R (equation (10)), the binding constant 
for repressor binding at the GS gene. Figure 4 shows the phase portrait corresponding to equations (8) 
and (9) with the parameter values f = 30, r] = 30, a = 0.002857, (5 = 0.002857 and \l = 0. The 
system is bistable for the parameter values quoted. The nullclines, obtained by putting ^ = 0,^ = 0, 
intersect at three points, the fixed points of the dynamics. The lower and upper fixed points are stable 
whereas the intermediate fixed point is unstable, in fact, a saddle node [30]. The stable manifold of the 
saddle node divides the uv— phase space into two basins of attraction. Trajectories starting in the lower 
(upper) basin of attraction end up at the lower (upper) stable fixed point as shown in figure 4. A trajectory 
initiated on the stable manifold stays on it and ends at the saddle node. A typical trajectory asymptotically 
approaches the unstable manifold as t — >• oo. A trajectory is obtained by plotting the values of u and v 
at different time points, determined by solving equations (8) and (9). The arrow direction on a trajectory 
denotes increasing time. 



Meaning/Defining Formula 

Rate constants for basal protein synthesis 
Rate constants for activated protein synthesis 
Binding constants of protein dimers S2 and 
P2 

Binding constants for the binding of protein 
dimers S2 and P2 at the genes GPand GS 
Rate constants for protein degradation 
Binding constant for repressor dimer binding 
at gene GS; K R thus denotes the strength of 
repression 

77 = 4 s -, ra tio of activated and basal rate 

' Jso 

constants for synthesis of S proteins 

£ = j^i ra ti° °f activated and basal rate 

constants for synthesis of P proteins 

/i = (^) 2 K S K R ; with J so , 75 and K s kept 

fixed, /x can be varied by changing K R thus 

providing a measure of repression strength 

a = (^fKp P K s 

(3=(h^yK S sKp 

abbreviations defined in equations (3) and (4) 
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Figure 4: Phase portrait described by equations (8) and (9). The dark and light solid lines represent the 
nullclines intersecting at three fixed points. The stable and unstable fixed points are denoted by solid 
and empty circles respectively. The stable manifold divides the phase space into two basins of attraction. 
Some typical trajectories are shown with arrow directions denoting increasing time. 

Figures 5(a) shows the plot of £ versus n exhibiting regions of monostability and bistability. The 
parameter values are the same as before with a — f3 — 0.002857 and /i — 0. The regions of bistability, 
enclosed within the red and black curves, correspond tons — 1 an d n s = 2 respectively. The difference 
in the locations of the two loops in the logarithmic plots clearly shows that the bistable region is of greater 
extent when the gene copy number is reduced from two to one. The region of bistability is decreased in 
extent when autorepression is taken into account (Figure 5(b) with /i = 0.005). Figure 6 shows the fx — (3 
plot with the regions of bistability falling within the red (ns = 1) and black (ns = 2) curves respectively. 
The value of /j, is changed by varying K R (equation 10) with /j, = 0.08 K R . 

A major advantage of combining a double positive feedback loop operating between two genes with 
autorepression of the expression of one of the genes lies in dosage compensation [16]. This relates to the 
fact that the fall in steady state protein levels, brought about by a reduction in the gene copy number, is 
less when autorepression is included, compared to the case when there is no autorepression. A measure 
of dosage compensation is provided by the quantity G, termed percentage gain, defined as 

GM = Xl ^- X ^ = °K l00 (11) 

Xidi = o) 

where x\ denotes the steady state concentration of S proteins when the copy number of the GS gene, 
n s , is one. The parameter // is a measure of the repression strength. G is calculated by keeping the 
mean level of S proteins to be the same in the two cases /i = and /i ^ when n s = 2. This is 
achieved by adjusting the binding constant K PP contained in the parameter a in equation (10). The 
other parameter values are rj = £ = 30.0, Jso = Jpo = 2.0, 75 = 7p = 1.0 and f3 = 0.002857. 
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Figure 5: Plots of £ versus rj showing regions of monostability and bistability when the parameter /i, a 
measure of the autorepression strength, is zero (a) and 0.005 (b). The regions of bistability are enclosed 
within the red and black curves with gene copy numbers n P = 2, n s = 1 and n P = 2, n s = 2 
respectively. 

Figure 7 shows the plot of G versus /j, for the parameter values mentioned. As /j, increases from zero, 
there is initially a sharp increase in the value of G followed by a slower growth which ultimately leads 
to a near- saturation of G values. The results obtained in the deterministic approach provide insight 
on the advantages of autorepression in the non-occurrence of the genetic disorder MODY. The normal 
functioning of pancreatic /3-cells requires the HNF — la and HNF — 4a protein levels to be high, i.e., 
the two-gene system should be in the ON state. The genesis of MODY lies in a substantial fraction of the 
/?-cells being in the OFF state. This is brought about by mutations in the hnf — la and hnf — 4a genes 
giving rise to a fall in the steady state protein levels. In terms of the two-gene model studied by us, the 
monostable high state, in which the levels of the P and S proteins are both high, represents the ON state 
of normal /3-cells. The system may enter a region of bistability, in which both the ON and OFF states are 
possible, due to the loss of a gene copy brought about by mutations. We will show in the next section 
that fluctuations in the protein levels are responsible for transitions between the ON and OFF states. In 
the deterministic scenario, the major advantages of the autorepressive feedback loop appear to be dosage 
compensation (figure 7) as well as a lesser possibility of the system being in the bistable region due to a 
reduction in gene copy number. The continuance of the system in the monostable high state ensures the 
normal functioning of cells. Similar conclusions are reached if the gene copy number n P is reduced from 
two to one. There is, however, an asymmetry in the S and P protein levels as the expression of the gene 
GP is not autorepressed. 
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Figure 6: /j, — (3 phase diagram with the regions of bistability falling within the red (n P — 2, ng — 1 ) 
and black (n P = 2, n s = 2) curves respectively. 
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Figure 7: Plot of percentage gain G (equation 11) versus /i, a measure of the autorepression strength. 

3 Stochastic Approach 

Consider the two-gene network to be originally in the monostable high state. In the deterministic for- 
malism, the system continues to be in the high, i.e., ON state even if it enters the region of bistability 
due to the loss of a gene copy. This is due to history dependence, since the system is initially in the 
ON state it continues to be in the ON state in the bistable region. The protein levels corresponding to 
the ON state are, however, lower in magnitude in the bistable region. In the pancreatic /3-ceils, the oc- 
currence of MODY is possible only when a sizable fraction of cells is in the OFF state. The ON^OFF 
and OFF^ON transitions can be understood only when stochasticity in gene expression is taken into 
account. We now give a simple physical picture of the origin of stochastic transitions [6]. In the case 
under consideration, the dynamical variables are the protein concentrations u and v. In the case of deter- 
ministic time evolution, trajectories starting in individual basins of attraction stay confined to the specific 
basins with no possibility of a trajectory crossing from one basin to another. In the stochastic approach, 
the trajectories are no longer deterministic as the dynamical variables u(t) and v(t) are fluctuating. In 
the deterministic case, given the initial state defined by (u(t = 0),v(t = 0)), the trajectory is fixed in the 
uv-phase space. In the stochastic case, different trajectories are generated in repeated trials. A transient 
fluctuation, if sufficiently strong, switches the system dynamics from one basin to the other brought about 
by the excursion of the trajectory across the boundary separating the two basins of attraction. In terms of 
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Figure 8: Distribution of steady state GS protein levels, P(u), in an ensemble of 4500 cells for repressor 
strengths (a) /i = 0.000027, (b) n = 0.00002, (c) n = 0.00001 and (d) fj, = 0.0 respectively. 

the pancreatic /3-cells, the switch to the OFF state hampers the normal functioning of the cells. 

For proper regulatory functions as transcription factors, the HNF — la and HNF — 4a protein 
levels are high with an optimal value as excessive protein amounts are known to be harmful rather than 
beneficial [16]. In this context, it is pertinent to undertake a comparison of the functional characteristics 
of two-gene network models with and without the autorepressive loop and with the mean protein levels 
kept at the same high values in the two cases. The last condition ensures the normal functioning of the 
cells in both the cases. In section 2, we have identified certain advantages of the autorepressive loop 
as regards the system dynamics in a deterministic framework. Our goal is now to identify the desirable 
features of the model incorporating both a double positive feedback loop and an autorepressive loop 
taking the stochastic aspects of the dynamics into consideration. This is done with the help of a detailed 
computer simulation based on the Gillespie algorithm [31]. The algorithm enables one to keep track 
of the stochastic time evolution of the system. The different biochemical reactions considered in the 
simulation are depicted in figures 2(a) and 2(b). The reactions are sixteen in number and are given by 

GS + PI -> GS* (12) 

GS* -> S (13) 
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GS -> S (14) 

5 + 5^52 (15) 

GS* ^GS + P2 (16) 

52^5 + 5 (17) 

S^<S> (18) 

CP + 52 -> GP* (19) 

GP* -> P (20) 

GP -> P (21) 

P + P -> P2 (22) 

GP* -> GP + 52 (23) 

P2 -> P + P (24) 

P -> $ (25) 

G5 + 52 -> G552 (26) 

G552 -> G5 + 52 (27) 



The different symbols are as explained in section 2. The stochastic rate constants, associated with the 
reactions, are C(i), i=l,...,16, in the appropriate units. The results of the simulation are shown in figures 
8-9. Figures 8(a)-(d) show the distribution of GS protein levels, P(u), in an ensemble of 4500 cells 
for repressor strengths \x = 0.000027, 0.00002, 0.00001 and 0.0 respectively after a simulation time of 
tmax = 2000 time units. The gene copy numbers are rip = 2 and n s = 1 so that the system is in the 
region of bistability. The values of the stochastic rate constants are G(2) = 56.0, G(3) = 2.0, G(4) = 4.0, 
G(5) = 280.0, G(6) = 100.0, G(7) = 1.0, G(8) = 10.0, G(9) = 50.0, G(10) = 2.0, G(ll) = 
4.0,G(12) = 280.0, G(13) = 100.0, C(U) = 1.0, G(15) = 10.0. The value of \i is changed by varying 
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the stochastic rate constant C(16). The value of the rate constant C(l) is changed to keep the mean 
protein levels to be the same when n s = 2, n P = 2 for all values of \i. The value ji — implies that 
only the double positive feedback loop contributes to the dynamics. The distribution P(u) is found to be 
bimodal, i.e., has two distinct peaks corresponding to the OFF and ON states. In all the cases, the cells 
are in the ON state at time t = 0. One finds that the fraction of cells in the OFF state decreases as the 
value of /i increases. In fact, when /i — 0, the number of cells which are in the OFF state is larger than 
that in the ON state. Since initially all the cells are in the ON state, a large number of ON — > OFF state 
transitions occur during the simulation time. For fi = 0, the reverse transition is, however, much rarer. 
The role of the autorepressive loop thus appears to be to reduce the number of stochastic transitions from 
the ON to the OFF state. This makes the occurrence of MODY, brought about by a sizeable fraction 
of the cell population existing in the OFF state, less probable. There are two distinct time scales over 
which protein fluctuations occur. The probability distribution P(u) versus u has a two-peaked structure. 
Fluctuations on a short time scale confine the u values to lie predominantly within individual peaks. 
The long time scale corresponds to the time at which large fluctuations occur bringing about transitions 
between states belonging to different peaks. The "escape time" is often very large and a quantitative 
measure is provided by the mean first passage time r [32]. In the present case, the values of t on _> off 
and toff-tON are quite large for different values of /i. The maximum simulation time tmax is 2000 time 
units for all values of /i. For fi = 0, t on ^ off is around 1000 time units whereas t off _, on is even larger. 
As fj, increases, t n^off increases whereas t ff^on decreases. For fi = 0.0005, t n->off is as large 
as 10 7 s. Because of large escape times, the probability distribution P(u) versus u is metastable on a large 
time scale [32]. Over shorter periods of time, the shape of the distribution remains more or less invariant. 

The plots in figure 9 are obtained for an ensemble of 4500 cells. For gene copy numbers n P = 2 
and n s = 2, the mean protein levels are adjusted to be the same irrespective of the values of . The 
parameter values are so chosen that the system is in the monostable high region. On reduction of ns to 
1 (one copy of the GS gene), the system enters the region of bistability and is in the ON state at time 
t = 0. After a period T = 2000 time units of stochastic time evolution, the percentage of cells in the 
OFF state is determined. The red curve shows this percentage as a function of the repression strength \x. 
The drop in the percentage of cells in the OFF state is found to be exponential. The black curve shows 
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Figure 9: For gene copy numbers rip = 2, ns = 1 and after a time interval T = 2000 time units 
of stochastic time evolution, the percentage of cells in an ensemble of 4500 cells in the OFF state (red 
curve) versus the repression strength ix with all the cells being in the ON state at time t — 0. The black 
curve shows the percentage of cells in the ON state versus ix with all the cells being in the OFF state at 
t = 0. 

the percentage of cells in the ON state after T = 2000 time units, with all the cells being initially in the 
OFF state. One finds that with increasing ix, the fraction of cells in the ON state becomes larger. The 
autorepressive loop has the effect of making the ON state more stable and the OFF state more unstable. 
This feature enhances the probability of the nonoccurrence of MODY as there are infrequent transitions 
from the ON to the OFF state. On the other hand, the system has a lesser probability of remaining stuck 
in the OFF state compared to the case when there is no autorepressive loop. 

4 Conclusion and Outlook 

In this paper, we have studied the functional characteristics of a motif consisting of a double positive 
feedback loop operating between two genes and a negative feedback loop in which the protein product 
of one gene represses its own synthesis. The motif appears in the gene regulatory network controlling 
the pancreatic /3-ceil function [16]. Loss of a gene copy due to mutations has been shown [16] to be 
responsible for abnormal /5-cell function resulting in MODY. We have studied the effect of reduced gene 
copy number on the dynamics of the model describing the two-gene motif. In a deterministic formalism 
based on differential rate equations, we identified regions of bistability in appropriate parameter regions. 
The stable steady states, designated as the OFF and ON states, correspond respectively to low and high 
protein levels. The normal /3-cells are expected to be in the monostable ON state. The occurrence of 
MODY is brought about by a fraction of /5-cells being in the OFF state. The ON^OFF switch can occur 
only in the bistable region. Negative feedback reduces the extent of the bistable region making it less 
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likely that the cellular state falls in this part of the phase diagram. The region of bistability, however, 
increases in size on reduction of the gene copy number making the ON^OFF transitions more probable. 
Negative feedback also produces a mechanism of dosage compensation (figure 7). The results hold true 
for a wide range of parameter values. Since switching to the OFF state is detrimental, one would have 
thought, from an evolutionary point of view, that two genes which are constitutively ON would be more 
appropriate. In reality, the genes hnf — la and hnf — Aa form a positive feedback loop. Cross-regulation 
between the two genes is established when the pancreatic cells receive signals to differentiate [16]. The 
positive feedback loop provides a stable mechanism of gene expression since the two genes reset each 
other's activity to the functional state under physiological perturbations. This serves to self-perpetuate 
the activity of the two genes and their targets in the pancreatic (3— cells. Normal functioning of these cells 
requires both the protein levels to be high. The system of two genes that are constitutively ON are less 
robust under physiological perturbations since there is no resetting mechanism by which both the genes 
are in the functional ON state. The theoretical suggestions of bistability due to the existence of a positive 
feedback loop [16, 18], backed up by the results of our mathematical model, should be tested in actual 
experiments. 

The ON^OFF switch is brought about by protein fluctuations the origin of which lies in stochastic 
gene expression. Our major finding is that negative feedback makes the ON^OFF transitions less prob- 
able and the OFF^ON transitions more probable. Thus the function of the negative feedback appears 
to be to protect the normal /9-cell function since the cell is more likely to be in the ON state in this case. 
The asymmetric response to fluctuations prevents switching off and facilitates switching on of the high 
expression state. In the deterministic scenario, one finds that the difference between the ON state and 
the unstable steady state protein level increases as the autorepression strength is increased whereas the 
difference between the unstable steady state and OFF state protein levels decreases on increasing the au- 
torepression strength. This may explain the asymmetry in the ON^OFF and OFF^ON switches when 
stochasticity is taken into account. For moderate strengths of autorepression, the system is locked in the 
ON state for extremely long times. In our simulations, we did not encounter ON^OFF switches for very 
long trajectories (~ 10 7 seconds) with /i = 0.0005. This translates into lifetimes measured in years and 
explains the delayed onset of the diabetic phenotype [16]. The phenotype generally appears after several 
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years indicating that the activation of the ON^OFF switch is rare. The average age at which MODY is 
manifest could thus be dictated by the probability that a sufficient number of /9-cells is locked in the OFF 
state. We have considered the simplest form of negative autoregulation in our two-gene model. There 
are recent suggestion that negative autoregulation of the HNF-4o; gene in the pancreatic /5-cells may be 
more complex [33]. Also, the number of transcription factor binding sites of the two genes is not known 
with certainty. Cooperative binding at multiple sites is expected to promote the stability of the gene ex- 
pression states. Our two-gene motif constitutes a minimal model which seeks to explains the desirable 
features of combining a double positive feedback loop with an autorepressive loop vis-a-vis the normal 
functional activity of /9-cells. The insight gained from the model study is expected to provide a basis for 
the investigation of more complex cases. 
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Figure Captions 

Figl. The two-gene network model. The protein products of the GP and GS genes activate each 
other's synthesis. There is also an autorepressive loop in which the proteins of the GS gene repress their 
own synthesis. 

Fig2. The reaction kinetic scheme of the two-gene model. The meanings of the symbols are explained 
in the text. 

Fig3. u versus n curves showing bistability and hysteresis. The solid (dashed) lines represent stable 
(unstable) steady states for gene copy numbers (a) n P = 2, n s = 2 and (b) n P — 2, ns — 1 . The 
parameter /i, a measure of the autorepression strength, is zero. 

Fig4. Phase portrait described by equations (8) and (9). The dark and light solid lines represent the 
nullclines intersecting at three fixed points. The stable and unstable fixed points are denoted by solid 
and empty circles respectively. The stable manifold divides the phase space into two basins of attraction. 
Some typical trajectories are shown with arrow directions denoting increasing time. 

Fig5. Plots of £ versus n showing regions of monostability and bistability when the parameter //, a 
measure of the autorepression strength, is zero (a) and 0.005 (b). The regions of bistability are enclosed 
within the red and black curves with gene copy numbers n P = 2, n s = 1 and n P = 2, n s = 2 
respectively. 

Fig6. fi — f3 phase diagram with the regions of bistability falling within the red (n P — 2, ns — 1 ) 
and black (n P = 2, ns = 2) curves respectively. 

Fig7. Plot of percentage gain G (equation 11) versus /i, a measure of the autorepression strength. 

Fig8. Distribution of steady state GS protein levels, P(u), in an ensemble of 4500 cells for repressor 
strengths (a) fj, = 0.000027, (b) ft = 0.00002, (c) fj, = 0.00001 and (d) fj, = 0.0 respectively. 

Fig9. For gene copy numbers n P = 2, ns = 1 and after a time interval T = 2000 time units of 
stochastic time evolution, the percentage of cells in an ensemble of 4500 cells in the OFF state (red 
curve) versus the repression strength ji with all the cells being in the ON state at time t — 0. The black 
curve shows the percentage of cells in the ON state versus \x with all the cells being in the OFF state at 
t = 0. 
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